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Abstract  Non-equilibrium  and  equilibrium  molecular 
dynamics  simulations  are  employed  to  study  the  mechan¬ 
ical  response  of  soda-lime  glass  (a  material  commonly  used 
in  transparent  armor  applications)  when  subjected  to  the 
loading  conditions  associated  with  the  generation  and 
propagation  of  planar  shock  waves.  Particular  attention  is 
given  to  the  identification  and  characterization  of  various 
(inelastic-deformation  and  energy-dissipation)  molecular- 
level  phenomena  and  processes  taking  place  at  the  shock 
front.  The  results  obtained  revealed  that  the  shock  loading 
causes  a  2-4%  (shock  strength-dependent)  density 
increase.  In  addition,  an  increase  in  the  average  coordina¬ 
tion  number  of  the  silicon  atoms  is  observed  along  with  the 
creation  of  smaller  Si-0  rings.  These  processes  are  asso¬ 
ciated  with  significant  energy  absorption  and  dissipation 
and  are  believed  to  control  the  blast/ballistic  impact  miti¬ 
gation  potential  of  soda-lime  glass.  This  study  was  also 
aimed  at  the  determination  (via  purely  computational 
means)  of  the  shock  Hugoniot  (i.e.,  a  set  of  axial  stress  vs. 
density/specific- volume  vs.  internal  energy  vs.  particle 
velocity  vs.  temperature)  material  states  obtained  in  soda- 
lime  glass  after  the  passage  of  a  shock  wave  of  a  given 
strength  and  on  the  comparison  of  the  computed  results 
with  their  experimental  counterparts.  The  availability  of  a 
shock  Hugoniot  is  critical  for  construction  of  a  high 
deformation-rate,  large-strain,  high  pressure  material 


M.  Grujicic  (£3)  •  B.  Pandurangan  •  W.  C.  Bell 
Department  of  Mechanical  Engineering, 

Clemson  University,  241  Engineering  Innovation  Building, 
Clemson,  SC  29634-0921,  USA 
e-mail:  gmica@exchange.clemson.edu 

B.  A.  Cheeseman  •  P.  Patel  •  G.  A.  Gazonas 

Army  Research  Laboratory — Survivability  Materials  Branch, 

Aberdeen,  Proving  Ground,  MD  21005-5069,  USA 


model  which  can  be  used  within  a  continuum-level  com¬ 
putational  analysis  to  capture  the  response  of  a  soda-lime 
glass-based  laminated  transparent  armor  structure  (e.g.,  a 
military  vehicle  windshield,  door  window,  etc.)  to  blast/ 
ballistic  impact  loading. 

Introduction 

While  novel  transparent  ceramics  and  transparent  poly¬ 
meric  materials  are  being  increasingly  used  in  various 
blast/ballistic-impact  resistant  vehicle  transparent  struc¬ 
tures  (e.g.,  windshields,  door  windows,  viewports,  etc.), 
ballistic  glass  remains  an  important  material  component  in 
these  structures  [1,2].  The  main  reasons  for  this  have  been 
discussed  in  great  detail  in  our  recent  work  [3-6]. 

Typically,  the  development  of  new  glass-based  trans¬ 
parent  impact  resistant  structures  involves  extensive  pro¬ 
totyping  and  laboratory/field  experimental  testing.  While 
these  experimental  efforts  are  critical  for  ultimate  valida¬ 
tion  of  the  transparent  impact  resistant  structures,  they  are 
generally  expensive,  time-consuming,  and  involve 
destructive  test  procedures.  In  addition,  it  is  well-estab¬ 
lished  (e.g.,  [3-6])  that  the  effectiveness  and  reliability  of 
the  computation-based  modeling  and  simulation  approa¬ 
ches  is  strongly  correlated  with  the  fidelity  of  the  associ¬ 
ated  material  models.  In  other  words,  it  is  critical  that  these 
models  realistically  describe  deformation/fracture  response 
of  the  subject  material  (ballistic  glass,  in  the  present  case) 
under  high-rate,  large-strain,  high-pressure  loading  condi¬ 
tions  encountered  during  blast/ballistic  impact.  Therefore, 
one  of  the  main  objectives  of  this  study  is  to  explore  the 
possibility  for  using  molecular-level  computational  simu¬ 
lations  to  improve  fidelity  of  the  models  representing  the 
material  under  consideration.  Specifically,  phenomena  and 
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processes  accompanying  the  formation  and  propagation  of 
shock- waves  within  soda-lime  glass  will  be  investigated  to 
derive  the  associated  functional  relations  between  different 
material  state  variables. 

A  comprehensive  literature  review  carried  out  as  part  of 
our  prior  work  [6]  revealed  that  the  mechanical  behavior  of 
glass  is  modeled  predominantly  using  three  distinct 
approaches:  (a)  molecular-modeling  methods  [7-10], 

(b)  continuum-material  approximations  [3-6,  11-14],  and 

(c)  models  based  on  explicit  crack  representation  [15,  16]. 
In  addition,  this  review  established  that  only  the  molecular- 
level  models  are  capable  of  identifying  nanometer  length 
scale  phenomena  and  the  associated  microstructure  evolu¬ 
tion  processes  accompanying  high-rate  deformation  of 
soda-lime  glass.  This  is  the  main  reason  that  molecular- 
level  methods  and  tools  are  employed  in  this  study. 

As  mentioned  above,  the  main  objective  of  this  study  is 
to  investigate  various  shock-wave-related  phenomena 
using  molecular-level  computational  methods. 

A  shock  wave  (or  simply  a  shock)  is  a  wave  which 
propagates  through  a  medium  at  a  speed  higher  than  the 
sound  speed  and  its  passage  causes  an  abrupt  and  discon¬ 
tinuous  change  in  the  material  state  variables  (e.g.,  pres¬ 
sure,  internal  energy,  density,  temperature,  and  particle 
velocity).  The  magnitude  of  the  state-variable  changes  and 
the  shock  speed  increases  with  the  strength  of  the  shock.  In 
reality,  since  no  material  can  support  a  discontinuity  in  its 
state  variables,  shocks  manifest  themselves  as  continuous 
(structured)  waves  with  a  steep  and  narrow  wave  front. 
While  acoustic  waves  give  rise  to  isentropic  changes  in  the 
material  state  variables,  passage  of  a  shock  is  typically 
associated  with  irreversible  (entropy-increasing)  changes  in 
the  same  variables.  The  reason  behind  this  difference  is 
that  shock  involves  very  high  strain  rates  that  bring  energy- 
dissipative  viscous  effects  into  prominence. 

Based  on  an  open-literature  review,  conducted  as  part  of 
this  study,  it  was  established  that  molecular-level  computa¬ 
tional  methods  were  first  employed  to  study  shocks  more 
than  30  years  ago  [17-20].  While  the  initial  focused  on 
shock  phenomena  in  dense  fluids,  subsequent  analyses  were 
extended  to  single-crystalline  solids.  The  key  findings  rela¬ 
ted  to  the  generation  and  propagation  of  shocks  in  the  single¬ 
crystalline  solids  investigated  can  be  summarized  as  follows: 
(a)  shock-wave  phenomena  are  inherently  more  complex  in 
solids  than  in  fluids,  because  solids,  in  addition  to  the  lattice 
parameter,  introduce  anew  length  scale  (i.e.,  the  size/spacing 
of  defects)  which  tends  to  control  the  nature/extent  of 
inelastic-deformation  processes  (mainly  responsible  for 
energy  dissipation  at  the  shock  front),  (b)  initial  (un-steady) 
shock  waves  tend  to  become  steady  (time-invariant)  as  they 
propagate  and  this  process  is  facilitated  by  the  transverse 
displacement  of  atoms  which  produce  inelastic  deformations. 


These  deformation  involves  concerted  slippage  of  atoms 
over  each  other  and  is  not  dominated  by  viscous  flow  as  in  the 
case  of  shocks  in  fluids,  and  (c)  in  order  to  eliminate  free- 
surface  effects  and  model  bulk-material  behavior,  molecu¬ 
lar-level  modeling  of  shock  is  typically  carried  out  using 
computational  systems  with  periodic  boundary  conditions 
(at  least  in  directions  transverse  to  the  shock-wave  propa¬ 
gation  direction).  As  a  consequence  of  the  use  of  the  periodic 
boundary  conditions  within  shock-wave  molecular-level 
simulations,  a  single  shock  (or  often  a  pair  of  converging 
shocks)  is  associated  with  each  computational  cell.  To  attain 
steady  wave  conditions  of  the  shock,  computational  domains 
sufficiently  long  in  the  direction  of  shock-wave  propagation 
have  to  be  employed.  For  a  typical  5000-7000  m/s  range  of 
the  shock  speed  and  a  typical  20  ps  computation  time  the 
minimum  computational-cell  length  required  is  of  the  order 
of  100-150  nm.  It  should  be  also  noted  that,  lateral  dimen¬ 
sions  of  the  computational  domain  have  to  be  also  suffi¬ 
ciently  large  to  avoid  spurious  effects  associated  with  the  use 
of  the  periodic  boundary  conditions.  Consequently,  com¬ 
putational  domains  involving  several  tens  of  thousands  of 
atoms  have  to  be  employed  in  shock-wave  simulations 
requiring  prolonged  wall-clock  simulation  times.  Further¬ 
more,  since  weak  and  moderate  shocks  in  single  crystalline 
solids  can  have  width  on  the  order  of  100-500  nm  they  could 
not  be  readily  modeled  using  the  molecular-level  methods. 
However,  as  will  be  shown  in  this  study,  when  the  energy 
dissipative  mechanisms  are  not  very  potent  (as  is  the  case  in 
soda-lime  glass),  moderate  shocks  are  associated  with  sig¬ 
nificantly  smaller  width  and  they  could  be  studied  using 
molecular-level  methods. 

As  mentioned  earlier,  a  preliminary  molecular-level 
computational  study  of  shock  generation/propagation  in 
soda-lime  glass  is  carried  out  in  this  study.  The  study  has 
the  following  two  main  objectives: 

(a)  To  determine  the  shock  Hugoniot  (centered  on  the 
initial  stress-free  quiescent  state)  of  soda-lime  glass. 
A  Hugoniot  is  the  locus  of  axial  stress  versus  specific 
volume  versus  energy  density  (vs.  particle  velocity  vs. 
shock  speed)  shocked  (“upstream”)  material  states. 
The  Hugoniot  is  often  used  in  the  derivation  of  the 
continuum-level  material  models  of  the  type  employed 
in  transient  non-linear  dynamics  analyses  of  the 
response  of  structures  to  shock  loading.  In  situations 
in  which  one  is  interested  only  in  the  problem  of  planar 
shock  propagation/interaction  (in  the  presence  of  uni¬ 
axial  strain  deformation  states),  a  complete  definition 
of  the  continuum-level  material  model  is  not  required. 
Instead,  the  knowledge  of  the  corresponding  Hugon- 
iots  (i.e.,  Hugoniot  functional  relations)  is  sufficient; 
and 


Springer 


7300 


J  Mater  Sci  (2011)  46:7298-7312 


(b)  The  Hugoniot  relations  mentioned  in  (a)  provide  a 
global  statement  of  mass,  momentum,  and  energy 
conservation  accompanying  shock-induced  material 
transition  from  a  given  initial  (“downstream”)  equi¬ 
librium  state  to  all  possible  final  (“upstream”) 
equilibrium  states  for  steady  planar  shock  waves  (of 
different  strengths).  However,  these  relations  provide 
no  information  about  the  structure  of  the  shock  front 
or  the  nature  of  the  dissipative  structural  rearrange¬ 
ment  mechanisms  that  lead  to  the  formation  of  a 
steady  shock  wave.  Hence,  the  second  objective  of 
this  study  is  to  carry  out  a  detailed  examination  of  the 
downstream,  shock-front,  and  upstream  material 
states  (as  represented  by  the  local  stresses,  strains, 
densities/specific  volumes,  temperatures,  etc.)  and 
molecular-level  morphology  to  identify  and  charac¬ 
terize  these  processes. 

The  organization  of  the  article  is  as  follows:  A  brief 
description  of  the  molecular-level  microstructure  of  glass 
including  its  random-network  representation  is  presented  in 
second  section.  A  brief  overview  of  the  molecular-level 
computational  procedure  including  the  computational  cell 
construction,  force  field  identification,  computational 
method(s)  selection,  shock-wave  generation,  and  the 
problem  definition  are,  respectively  presented  in  third 
section.  The  key  results  obtained  in  this  study  are  presented 
and  discussed  in  fourth  section.  A  summary  of  the  study 
carried  out  and  the  key  results/conclusions  are  given  in  fifth 
section. 


Molecular-level  microstructure  of  glass 

One  of  the  basic  properties  of  all  glass-like  materials  is 
their  lack  of  long-range  atomic  order  which  classifies 
them  as  amorphous  materials.  For  instance,  the  atomic 
arrangement  in  pure  silicate  glass  (i.e.,  fused  silica)  is 
highly  random  relative  to  the  crystalline  Si-0  modifica¬ 
tions  like  quartz  or  cristobalite.  To  describe  the  random 
atomic  arrangement  within  glass  the  so-called  random 
network  model  [21]  is  typically  employed.  Such  a  model 
represents  amorphous  materials  as  a  three-dimensional 
linked  network  of  polyhedra  with  central  cations  of 
various  coordinations  depending  on  the  character  of  the 
atomic  constituents.  In  the  case  of  silicate-based  glasses 
like  fused  silica  and  soda-lime  glass,  the  polyhedral 
central  cation  is  silicon.  In  this  case,  each  silicon  is 
surrounded  by  four  oxygen  anions  and  forms  a  Si044- 
tetrahedron,  whereby  each  oxygen  is  bonded  to  (or 
bridges)  two  silicon  atoms.  Other  polyhedra  may  exist  in 
the  network  depending  on  the  valence  of  the  added 
polyhedral-center  cations. 


Ceramic  glasses  like  soda-lime  glass  investigated  in  this 
study  are  generally  a  mixture  of  various  oxides.  Oxide 
constituents  of  the  ceramic  glasses  are  generally  classified 
as:  (a)  Network  formers  oxides  (e.g.,  Si02,  B203,  Ge02) 
falling  into  this  category  have  the  tendency  to  form  a 
continuous  network  with  “bridging”  oxygen  atoms;  and 
(b)  Network  modifiers  oxides  falling  into  this  category  are 
typically  based  on  alkali  (e.g.,  Na,  K,  etc.)  and  alkaline 
earth  (e.g.,  Ca,  Mg,  etc.)  metals.  In  the  specific  case  of 
soda-lime  glass,  investigated  in  this  study,  the  14  wt% 
Na20  and  9  wt%  CaO  additions  pure  Si02  glass  both  act  as 
network  modifiers.  When  alkali/alkaline  earth  metal -based 
oxides  are  added  to  pure  Si02  glass,  excess  oxygen  which 
results  from  the  ionic  dissociation  of  these  oxides  is 
incorporated  into  the  glass  network.  This  process  disrupts 
network  continuity  since  it  replaces  single  continuity¬ 
forming  bridging  oxygen  atom  with  two  (network-break¬ 
ing)  singly  charged  non-bridging  oxygen  ions.  The  latter 
oxygen  ions  are  the  result  of  charge  transfer  from  the  ini¬ 
tially  double  charged  oxygen  ions  (produced  during  dis¬ 
sociation  of  alkali/alkaline  earth  metal-based  oxides)  to  the 
initially  neutral  bridging  oxygen  atoms.  The  metallic  cat¬ 
ions  formed  during  dissociation  tend  to  hover  around  the 
non-bridging  oxygen  ions  for  local  charge  neutrality. 

Within  the  random  network  model,  the  micro  structure  of 
glass  is  described  using  several  network  state  parameters. 
Among  these,  the  most  frequently  used  are:  (a)  the  so- 
called  R  parameter,  defined  as  the  average  number  of 
oxygen  ions  per  network  forming  cation.  In  a  glass  of  a 
given  chemistry,  this  parameter  is  effectively  micro¬ 
structure-insensitive  since  its  value  depends  solely  on  the 
chemical  composition  of  glass.  For  example,  in  the  case  of 
fused  silica,  in  which  there  are  two  oxygen  anions  for  every 
network  forming  silicon  cation,  the  R  value  is  2.0.  On  the 
other  hand,  in  the  case  of  soda-lime  glass,  the  introduction 
of  additional  oxygen  ions  (without  the  introduction  of 
additional  network  forming  cations)  causes  the  R  value  to 
increase  to  ca.  2.41.  In  accordance  with  the  earlier  dis¬ 
cussion  regarding  the  role  of  glass  modifiers,  it  is  apparent 
that  a  larger  value  of  the  R  parameter  implies  a  more  open, 
weaker  glass  network.  The  effect  of  network  formers  on  the 
R  parameter  is  more  complicated  and  depends  on  the  net¬ 
work  former  coordination  number  as  well  as  its  concen¬ 
tration.  Overall  the  R  parameter  is  a  measure  of  the 
maximum  extent  of  connectivity  that  can  be  achieved 
within  a  glass  of  a  given  chemistry.  It  should  be  noted  that 
since  shock  loading  is  assumed  to  only  produce  physical/ 
microstructural  changes  in  glass,  the  R  value  is  not 
expected  to  change  during  shock  loading;  (b)  the  so-called 
X  parameter  which  defines  the  average  number  of  non¬ 
bridging  (connected  to  only  a  single  network  forming 
cation)  oxygen  atoms  per  network  polyhedron.  The  equi¬ 
librium  value  of  this  parameter  can  again  be  calculated 
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using  the  known  chemical  composition  of  glass.  For 
example,  this  parameter  takes  on  a  zero  equilibrium  value 
in  the  case  of  fused  silica  since  all  oxygen  ions  are  of  a 
bridging  type.  However,  in  the  case  of  soda-lime  glass  for 
the  reasons  mentioned  above,  the  equilibrium  value  of  this 
parameter  is  ca.  0.81.  It  should  be  noted  that  shock  loading- 
induced  damage  in  glass  will  generally  increase  the  number 
of  non-bridging  oxygen  ions  (through  Si-0  bend  breaking), 
resulting  in  an  increase  of  the  X  parameter  beyond  its 
equilibrium  value;  and  (c)  the  so-called  Y  parameter  which 
defines  the  average  number  of  bridging  (connected  to  two 
network  forming  cations)  oxygen  atoms  per  network 
polyhedron.  The  equilibrium  values  for  this  parameter  in 
fused  silica  and  soda-lime  glass  are  4.0  and  3.19,  respec¬ 
tively.  Shock  loading  will  generally  lower  the  value  of  this 
parameter  below  the  equilibrium  value. 

In  addition  to  chemical  modifications  of  glass,  changes 
in  the  micro  structure  of  this  material  can  be  brought  about 
by  mechanical  loading/deformation  (typically  requiring 
several  GPa  pressure  levels).  Specifically,  high  pressure 
may  result  in  a  reorganization  of  the  atomic  network  (phase 
change)  in  the  form  of  changes  to  the  coordination  of  the 
network  forming  cations.  These  phase  changes  can  be  of 
first  order,  which  are  characterized  by  the  formation  of  a 
distinct  high-pressure  phase  at  a  nominally  constant  pres¬ 
sure,  or  they  may  be  of  second  order  which  are  phase 
changes  which  involve  a  continuous  morphing  of  the 
original  phase  into  the  final  high-pressure  phase  over  a 
range  of  pressures.  These  phase  transformations  may  be 
associated  with  significant  volume  changes  and,  since 
phase-transformation  induced  energy  absorption  is  a  well- 
documented  phenomenon  responsible  for  high  toughness 
levels  in  TRIP  steels  and  partially  stabilized  crystalline 
ceramics,  it  is  of  interest  to  this  study. 

The  phase  transformations  of  glass  which  are  investi¬ 
gated  in  this  study  occur  in  the  range  of  3-5  GPa.  The  phase 
transformations  in  question  have  not  been  fully  investigated 
or  characterized  in  terms  of  the  random  network  model  or 
using  other  micro  structural  parameters.  However,  the 
presence  of  these  transformations  has  been  inferred  indi¬ 
rectly  through  the  changes  in  the  shock  Hugoniot  relations 
for  soda-lime  glass  [22].  These  changes  typically  result  in 
an  anomalous  stiffness-decreasing  response  of  the  material 
during  loading.  This  pressure  range  was  chosen  for  the 
current  investigation  as  it  is  consistent  with  the  levels  seen 
in  glass  structures  subject  to  ballistic/blast  loading  and  is 
associated  with  relatively  modest  (3-7%)  volume  changes. 
These  phase-transformations  should  not  be  confused  with 
those  seen  to  take  place  at  substantially  higher  pressures 
(ca.  >20  GPa),  which  are  considered  to  be  of  second  order 
and  associated  with  substantially  larger  volume  reductions 
and  with  the  formation  of  stishovite,  an  octahedrally  coor¬ 
dinated  glass  phase. 


Molecular-level  analysis  of  soda-lime  glass 

As  mentioned  earlier,  molecular-level  computational 
methods  have  been  employed  in  this  study  to  investigate 
shock  formation  and  propagation  in  soda-lime  glass. 
Within  these  methods,  all  atoms,  ions,  and  bonds  are 
explicitly  accounted  for  and  molecular  mechanics, 
dynamics  or  Monte  Carlo  algorithms  are  used  to  quantify 
the  behavior  of  the  material  under  investigation. 

While  ab  initio  quantum  mechanics  methods  have  the 
advantage  over  the  molecular-level  methods  since  they  do 
not  require  parameterization,  they  have  a  serious  short¬ 
coming.  Namely,  due  to  prohibitively  high  computational 
cost,  they  can  be  currently  employed  only  for  systems 
containing  no  more  than  a  few  hundred  atoms/particles.  As 
will  be  shown  below,  while  ab  initio  quantum-mechanics 
calculations  are  not  directly  used  in  this  study,  some  of  the 
computational  ab  initio  quantum  mechanics  results  are 
used  in  the  parameterization  of  the  material  model  at  the 
molecular  length/time  scale.  Utility  of  the  molecular-level 
computational  results  is  greatly  dependent  on  accuracy  and 
fidelity  of  the  employed  force  field  (the  mathematical 
expressions  which  describe  various  bonding  and  non¬ 
bonding  interaction  forces  between  the  constituents  of  the 
molecular-level  model).  In  this  study,  the  so-called  Con¬ 
densed-phase  Optimized  Molecular  Potentials  for  Atomis¬ 
tic  Simulation  Studies  (COMPASS)  force  field  is  used 
[23,  24].  This  highly  accurate  force  field  is  of  an  ab  initio 
type  since  most  of  its  parameters  were  determined  by 
matching  the  predictions  made  by  the  ab  initio  quantum 
mechanics  calculations  to  the  condensed-matter  experi¬ 
mental  data.  Hence,  it  should  be  recognized  that  the 
COMPASS  force  field  is  a  prime  example  of  how  the 
highly  accurate  results  obtained  on  one  length/time  scale 
(quantum  mechanic/electronic,  in  the  present  case)  and  the 
experimental  data  can  be  combined  to  parameterize  mate¬ 
rial  models  used  at  coarser  length/time  scale  (the  molecular 
length/time  scale,  in  the  present  case). 

Formulation  of  a  molecular-level  simulation  problem 
requires,  at  a  minimum,  specification  of  the  following  three 
aspects:  (a)  a  molecular-level  computational  model  con¬ 
sisting  of  atoms,  ions,  functional  groups,  and/or  molecules; 
(b)  a  set  of  force  field  functions;  and  (c)  a  computational 
method(s)  to  be  used  in  the  simulation.  More  details  of 
these  three  aspects  of  the  molecular-level  modeling  and 
simulation  of  soda-lime  glass  are  provided  below. 

Computational  model 

At  the  molecular  level,  soda-lime  glass  is  modeled  as  a 
discrete  material  consisting  of  (a)  silicon  (Si)  and  oxygen 
(O)  atoms  mutually  bonded  via  a  single  covalent  bond  and 
forming  a  connected,  non- structured/amorphous  network  of 
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silica  (Si044-)  tetrahedra;  (b)  oxygen  anions  (O2-) 
attached  as  terminal  functional-groups  to  the  fragmented 
silica  tetrahedra  network;  and  (c)  sodium  cations  (Na+) 
dispersed  between  fragmented  silica  tetrahedra  networks 
and  ionically  bonded  to  the  oxygen  anions. 

While  glass  is  an  amorphous  material  and  does  not 
possess  any  long-range  regularity  in  its  atomic/molecular 
structure,  modeling  of  bulk  behavior  of  glass  is  typically 
done  at  the  molecular  level  by  assuming  the  existence  of  a 
larger  (amorphous)  unit  cell.  Repetition  of  this  cell  in  the 
three  orthogonal  directions  (the  process  also  known  as 
application  of  the  “periodic  boundary  conditions”)  results 
in  the  formation  of  an  infinitely  large  bulk-type  material. 
This  procedure  was  adopted  in  this  study. 

The  parallelopiped- shaped  computational  cell  used  in 
this  study  contained  2,304  particles  (Si-512,  Na+-512, 
0— 256,  and  0-1024).  The  computational  cell  edge-lengths 
were  initially  set  to  a  =  1  All  nm  and  b  =  c  =  1.868  nm, 
yielding  a  soda-lime  glass  initial  nominal  density  of 
2.613  g/ cm3.  The  three  edges  (< a ,  b ,  and  c)  of  the  cell  were 
aligned,  respectively,  with  the  three  coordinate  axes  (x,  y, 
and  z). 

To  create  the  initial  particle  configuration  in  the  unit 
cell,  the  Visualizer  [25]  program  from  Accelrys  was  first 
used  to  construct  a  short  silica-chain  fragment.  The  frag¬ 
ment  was  then  “grown”  by  a  duplicate-and-attach  process 
using  the  same  program.  The  resulting  silica  network 
(along  with  additional  sodium  cations  and  oxygen  anions) 
was  next  used  within  the  Amorphous  Cell  program  [26] 
from  Accelrys  to  randomly  populate  the  computational  cell 
while  ensuring  that  the  target  material  density  of  2.613  g/ 
cm3  was  attained.  An  example  of  a  typical  molecular-level 
topology  within  a  single  unit  cell  is  displayed  in  Fig.  1. 

Force-fields 

As  stated  above,  the  behavior  of  a  material  system  at  the 
molecular-level  is  governed  by  the  appropriate  force-fields 
which  describe,  in  an  approximate  manner,  the  potential 
energy  hyper- surface  on  which  the  atomic  nuclei  move.  In 
other  words,  the  knowledge  of  force-fields  enables  deter¬ 
mination  of  the  potential  energy  of  a  system  in  a  given 
configuration.  In  general,  the  potential  energy  of  a  system 
of  interacting  particles  can  be  expressed  as  a  sum  of  the 
valence  (or  bond),  cross-term  and  non-bond  interaction 
energies.  Each  of  these  energy  components,  in  turn,  con¬ 
tains  a  number  of  subcomponents.  For  example,  the 
valence  energy  contains  the  contributions  from  single  bond 
stretching,  two-bond  angle  changes,  three-bond  dihedral 
bond  angle  changes,  etc.  The  cross  energy  term  includes 
the  contributions  from  the  interactions  of  the  subcompo¬ 
nents  of  the  valence  terms  (e.g.,  stretch/stretch  interaction, 
stretch/bond-angle  change  interaction,  etc.).  The  non-bond 


Fig.  1  a  The  computational  unit  cell  for  soda-lime  glass  molecular- 
level  simulations  used  in  this  study;  and  b  an  example  of  the  local 
atomic  structure  (Color  figure  online) 

interaction  term  accounts  for  the  interactions  between  non- 
bonded  atoms  and  includes  the  van  der  Waals  energy  and 
the  Coulomb  electrostatic  energy. 

As  mentioned  earlier,  the  present  molecular-level  anal¬ 
ysis  of  soda-lime  glass  employs  the  COMPASS  [23,  24] 
force-field  for  various  bond  and  non-bond  interaction 
energies.  A  summary  of  the  COMPASS  force-field  func¬ 
tions  can  be  found  in  our  previous  study  [27]. 

Computational  method 

Both  molecular  statics  and  molecular  dynamics  simulations 
were  employed  in  this  study.  Within  the  molecular  statics 
approach,  the  unit-cell  potential  energy  is  minimized  with 
respect  to  the  position  of  the  constituent  particles/atoms. 
The  potential  energy  minimization  within  Discover  [28] 
(the  atomic  simulation  program  from  Accelrys  used  in  this 
study)  is  carried  out  by  combining  the  Steepest  Descent, 
Conjugate  Gradient,  and  the  Newton’s  minimization 
algorithms.  These  algorithms  are  automatically  inactivated/ 
activated  as  the  atomic  configuration  is  approaching  its 
energy  minimum  (i.e.,  the  Steepest  Descent  method  is 
activated  at  the  beginning  of  the  energy  minimization 
procedure  while  the  Newton’s  method  is  utilized  in  the 
final  stages  of  this  procedure). 

As  will  be  discussed  in  greater  detail  in  section  “Results 
and  discussion,”  molecular  statics  is  employed  to  deter¬ 
mine  the  state  of  the  material  swept  by  a  shock.  As  will  be 
shown,  this  procedure  is  based  on  the  use  of  (bonding  and 
non-bonding)  potential  energy  components  and  neglects 
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shock-induced  changes  in  the  (configurational)  entropy  of 
the  system.  To  assess  the  consequence  of  this  simplifica¬ 
tion,  the  approach  described  in  [29]  was  considered.  This 
approach  defines  a  dimensionless  parameter  and  states  that 
when  this  parameter  is  significantly  smaller  than  unity, 
entropy  effects  can  be  neglected.  Unfortunately,  detailed 
temperature  and  pressure  dependencies  of  the  material 
mechanical  response  of  soda-lime  glass  needed  to  evaluate 
this  parameter  were  not  available  with  sufficient  fidelity. 
Hence,  the  results  obtained  by  the  application  of  this  pro¬ 
cedure,  which  suggest  that  the  entropy  effects  may  not  be 
highly  critical,  cannot  be  accepted  with  a  high  level  of 
confidence. 

Within  the  molecular  dynamics  approach,  gradient  of 
the  potential  energy  with  respect  to  the  particle  positions  is 
first  used  to  generate  forces  acting  on  the  particles  and, 
then,  the  associated  Newton’s  equations  of  motion  (for  all 
particles)  are  integrated  numerically  to  track  the  temporal 
evolution  of  the  particle  positions.  Both  the  equilibrium 
and  the  non-equilibrium  molecular  dynamics  methods  are 
employed  in  this  study.  Within  the  equilibrium  molecular- 
dynamics  methods,  the  system  under  consideration  is 
coupled  to  an  (external)  environment  (e.g.,  constant  pres¬ 
sure  piston,  constant  temperature  reservoir,  etc.)  which 
ensures  that  the  system  remains  in  equilibrium  (i.e.,  the 
system  is  not  subjected  to  any  thermodynamic  fluxes).  As 
will  be  discussed  in  next  section,  NVT  [where  N  is  the 
(fixed)  number  of  atoms,  V,  the  computational  cell  volume 
(also  fixed),  and  T(= 298  K)  is  the  temperature]  equilibrium 
molecular  dynamics  is  employed  in  the  first  stage  of  the 
shock  generation  procedure.  In  addition,  as  will  be  dis¬ 
cussed  in  section  “Results  and  discussion”,  NVE  ( E  is  the 
total  energy)  equilibrium  molecular  dynamics  is  also 
employed  during  determination  of  the  shock  Hugoniot. 
Within  non-equilibrium  molecular  dynamics,  the  system  is 
subjected  to  large  perturbations  (finite  changes  in  the  axial 
parameter  of  the  computational  cell,  in  the  present  case) 
which  create  a  thermodynamic  flux  (i.e.,  the  flux  of  energy 
and  momentum,  in  the  present  case).  More  details  regard¬ 
ing  the  use  of  Discover  to  carry  out  molecular  statics  and 
molecular  dynamics  analyses  can  be  found  in  our  prior 
study  [6]. 

Shock-wave  generation 

To  generate  a  planar  shock  (or  more  precisely  a  pair  of 
planar  shocks)  within  the  computational  cell,  the  following 
procedure  is  employed: 

(a)  At  the  beginning  of  the  analysis,  a  “sufficiently  long” 
NVT  molecular  dynamics  simulation  is  carried  out  to 
equilibrate  the  system/material. 


(b)  The  shock  is  then  initiated  (and  driven)  by  continu¬ 
ously  contracting  the  computational  cell  v-direction 
lattice  parameter  a  as: 

a(t)  =  a(t  =  0)  —  2  uvt  (1) 

where  t  denotes  time,  uv  is  the  so-called  piston  velocity 
(or  equivalently  the  particles  upstream  velocity)  in  the 
v-direction.  uv  is  varied  over  a  range  between  187.5  and 
1500  m/s  to  simulate  the  generation  and  propagation  of 
shock  of  various  strengths.  Meanwhile,  computational-cell 
transverse  lattice  parameters  b  and  c  are  kept  constant  to 
obtain  planar  (uniaxial-strain)  shock  conditions.  In  this 
process,  the  computational  cell  faces  normal  to  the  shock 
propagation-direction  behave  very  similarly  to  the  impact- 
surface  of  a  plate-like  target  subjected  to  a  so-called 
symmetric  “flyer-plate”  impact  test  [30].  The  procedure 
employed  here  generates  a  pair  of  shock  waves  which 
propagate,  at  a  shock  speed  Us,  from  the  cell  boundaries 
toward  its  center.  It  should  be  noted  that  while  the  proce¬ 
dure  employed  involved  contraction  of  the  unit-cell  at  only 
one  of  its  faces,  due  to  the  use  of  the  periodic  boundary 
conditions,  a  pair  of  shocks  was  generated.  In  other  words, 
it  was  not  possible  to  generate  a  single  shock  while 
maintaining  material  continuity  (i.e.,  satisfying  the  periodic 
boundary  condition).  As  schematically  shown  in  Fig.  2, 
these  shock  waves  leave  behind  a  “shocked”  material  state 
characterized  by  a  higher  material  density  (as  well  as 
internal  energy,  temperature,  stress,  particle  velocity,  and 
entropy). 

The  aforementioned  procedure  for  shock-wave  genera¬ 
tion  and  the  subsequent  molecular  statics/dynamics  analy¬ 
ses  are  carried  out  through  the  use  of  a  Discover  input  file 
[28]  which  is  written  in  a  Basic  Tool  Command  Language 
(BTCL).  This  enabled  the  use  of  a  scripting  engine  that 
provides  very  precise  control  of  simulation  jobs,  e.g.,  a  cell 
deformation  to  be  carried  out  in  small  steps  each  followed 
by  a  combined  energy-minimization/molecular-dynamics 
simulation  run. 

Problem  formulation 

The  problem  addressed  in  this  study  involved  generation  of 
shock  waves  of  different  strengths  (using  the  aforemen¬ 
tioned  computational  cell  parameter  contraction  method), 
determination  of  the  associated  shock-Hugoniot  relations 
and  identification  and  elucidation  of  the  main  molecular- 
level  inelastic-deformation/energy  dissipation  processes 
taking  place  at  or  in  the  vicinity  of  the  shock  front.  The 
procedure  for  shock-wave  generation  was  presented  in  the 
previous  section. 
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Fig.  2  A  schematic  of  the 
generation  of  a  pair  of  shocks  in 
a  molecular-level  system  via  the 
process  of  computational-cell 
parameter  contraction  (Color 
figure  online) 


As  far  as  the  shock  Hugoniot  determination  is  con¬ 
cerned,  it  entailed  the  knowledge  of  the  shock-wave  pro¬ 
files  (and  their  temporal  evolution)  for  the  axial  stress, 
material  density,  particle  velocity,  internal  energy,  and 
temperature.  The  latter  are  obtained  by  lumping  particles/ 
atoms  and  their  (bond  and  non-bond)  potential  and  kinetic 
energy  contribution,  into  fixed- width  bins,  in  the  order  of 
their  axial  coordinates.  As  will  be  shown  in  the  next  sec¬ 
tion,  two  types  of  bins  are  used:  (a)  a  Lagrangian-type 
which  is  fixed  to  the  initial/reference  state  of  the  compu¬ 
tational  cell  and  (b)  a  moving-type  which  is  attached  to  the 
advancing  shock  front. 

Identification  of  the  molecular-level  inelastic-deforma- 
tion/energy  dissipation  processes  entailed  a  close  exami¬ 
nation  of  the  changes  in  a  material  bond  structure  and 
topology  caused  by  the  passage  of  the  shock. 

Results  and  discussion 

Shock-wave  observation  and  structure  characterization 

An  example  of  the  typical  results,  obtained  in  this  study, 
pertaining  to  material  molecular-level  microstructure/ 
topology  evolution  caused  by  a  continuous  axial  contrac¬ 
tion  of  the  computational  cell  is  displayed  in  Fig.  3a-d. 
The  results  displayed  in  these  figures  clearly  reveal  the 


generation  of  a  pair  of  planar  shock  waves  at  the  two 
opposing  y-z  faces  of  the  computational  cell,  Fig.  3a,  and 
their  subsequent  propagation  toward  the  center  of  the  cell, 
Fig.  3b-d.  An  approximate  location  of  the  center-point  of 
the  two  shocks  is  indicated  using  arrowheads  in  Fig.  3a-d. 
The  results  displayed  in  Fig.  3a-d  show  that  the  shock 
waves  remain  fairly  planar  during  their  motion.  The  anal¬ 
ysis  of  shock-wave  propagation  was  terminated  at  a  time 
which  is  shorter  than  the  two- shock  collision  time.  How¬ 
ever,  the  analysis  can  be  extended  throughout  the  shock- 
interaction,  which  will  be  done  in  our  future  study. 

While  the  results  displayed  in  Fig.  3a-d  provide  clear 
evidence  for  the  formation  and  propagation  of  a  pair  of 
opposing  planar  shock  waves,  they  do  not  offer  any 
information  about  the  structure/shape  of  the  shock-wave 
front  or  about  the  state  of  the  (upstream)  material  swept  by 
the  shock.  The  latter  aspects  of  shock- wave  generation  and 
propagation  within  soda-lime  glass  are  addressed  in  the 
remainder  of  this  section  and  in  the  subsequent  sections. 

To  reveal  the  structure/shape  of  the  shock  wave,  the 
method  of  (Lagrangian)  bins  described  in  section  “shock- 
wave  generation”  is  employed.  In  this  case,  the  bins  are 
stationary  since  they  are  fixed  to  in  the  (initial)  reference 
atomic  configuration  of  the  computational  cell.  Conse¬ 
quently,  the  same  atoms  are  associated  with  a  given  bin 
throughout  the  entire  molecular  dynamics  simulation. 
These  bins  (with  an  increased  thickness,  for  clarity)  are 


Springer 


J  Mater  Sci  (2011)  46:7298-7312 


7305 


Fig.  3  Temporal  evolution  of 
the  molecular-level  material 
microstructure  accompanying 
generation  and  propagation  of  a 
pair  of  planar  shocks  in  soda- 
lime  glass  (Color  figure  online) 


depicted  in  Fig.  3  top  part,  the  initial  configuration. 
Examples  of  the  typical  results  obtained  through  the  use  of 
this  method  are  displayed  in  Fig.  4a,  b.  The  results  dis¬ 
played  in  these  figures  are  obtained  under  identical  con¬ 
ditions  except  for  the  rate  of  axial  contraction  of  the 
computational  cell  which  is  100%  higher  in  the  case  of 
Fig.  4b  (750  m/s)  relative  to  Fig.  4a  (375  m/s).  In  these 
figures,  particle  velocities  at  different  simulations  (i.e., 
post- shock- wave  generation)  times,  1  =  60  fs,  2  =  120  fs, 
3  =  180  fs,  etc.  are  plotted  against  the  Lagrangian  bin 
center  v-location  (for  the  shock  propagating  to  the  right, 
only).  Brief  examination  of  the  results  displayed  in  Fig.  4a, 
b  reveals  that: 

(a)  two  shock  waves  are  generated  (only  the  right- 
propagating  shock  is  shown  though)  at  the  computa¬ 
tional  cell  faces  normal  to  the  v-direction.  These 


shocks  subsequently  propagate  toward  the  computa¬ 
tional-cell  center; 

(b)  after  a  brief  transient  period,  the  shocks  appear  to 
attain  steady  wave  profiles  (i.e.,  a  nearly  time- 
invariant  profile  within  a  reference  frame  which  is 
attached  to,  and  moves  with,  the  shock  front).  This  is 
demonstrated  for  the  weaker  shock  case  in  Fig.  4c  in 
which  the  shock  profile  for  the  curves  labeled  3  and  4 
in  Fig.  4a  are  redrawn  within  a  moving  frame 
attached  to  the  inflection  point  of  the  shock  profile; 
and 

(c)  both  the  particle  velocity  and  the  shock  speed  increase 
with  the  computational-cell  contraction  rate.  It  should 
be  noted  that  the  curves  bearing  the  same  numerical 
label  in  Fig.  4a,  b  correspond  to  the  same  simulation 
time. 
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Fig.  4  Temporal  evolution  of 
the  particle  velocity  associated 
with  the  propagation  of  two 
approaching  shock  waves  under 
the  imposed  computational  cell 
contraction  rate  of  a  375  m/s; 
and  b  750  m/s.  The  simulation 
time  associated  with  each  of  the 
curves  is  equal  to  the  curve 
number  label  multiplied  by 
60  fs;  c  superposition  of  the 
curves  3  and  4  from  part 
a  suggesting  steady  nature  of 
the  shock  wave;  and 
d  superposition  of  the  curves  3, 
4,  and  5  from  parts  a  and 
b  demonstrating  that  the  shock 
width  decreases  while  shock 
speed  increases  with  an  increase 
in  shock  strength.  Note  that  the 
particle  velocities  are 
normalized  with  respect  to  the 
respective  maximum  value 
(Color  figure  online) 
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It  should  also  be  noted  that  no  thermostat  was  used  in  the 
present  non-equilibrium  molecular  dynamics  simulations, 
so  that  the  steady-wave  shock  profile  is  a  natural  conse¬ 
quence  of  a  balance  between  the  continuous  supply  of 
momentum  to  the  system  (through  the  continuous  compu¬ 
tational  cell  axial  contraction)  and  the  observed  lateral 
motion  of  the  atoms  in  the  continuously  enlarged  upstream 
material  domain  swept  by  the  shock.  It  should  be  noted  that, 
due  to  the  use  of  periodic  boundary  conditions  in  the  lateral 
directions,  despite  the  fact  that  the  unit  cell  was  constrained 
in  these  directions,  atoms  are  free  to  pass  through  the 
computational  cell  lateral  boundaries.  In  general,  the  use  of 
a  thermostat  modifies  the  (F  =  ma  Newtonian-type)  equa¬ 
tions  of  motion  solved  within  the  molecular  dynamics 
simulations  by  the  introduction  of  a  (velocity-proportional) 
viscous-dissipation  term.  It  is  well-established  that,  when 
shock  formation  and  propagation  is  analyzed  within  a  con¬ 
tinuum  framework,  the  use  of  a  viscous-dissipation  term  is 
mandatory  for  the  attainment  of  a  steady- wave  shock  pro¬ 
file.  This  fact  has  often  been  used  as  a  justification  for  the  use 
of  a  (local  or  global)  thermostat  within  molecular-level 
simulations  of  shock-wave  formation/propagation.  While 
such  practices  greatly  facilitate  the  attainment  of  a  steady 
shock,  they  cannot  be  readily  defended  since  shock 


formation  and  propagation  is  generally  considered  to  be  an 
adiabatic  (no  system/surrounding  energy  exchange)  process 
due  to  the  near-instantaneous  material  transition  to  the 
shocked  state.  In  addition,  shock  generation/propagation  is  a 
non-isentropic  process  due  to  the  attendance  of  various 
energy  dissipation  mechanisms.  These  were  the  determining 
factors  in  the  decision  to  not  use  a  thermostat  in  this  study. 
Due  to  the  lack  of  a  thermostat,  the  molecular  dynamics 
simulations  employed  in  this  study  can  be  characterized  as 
being  of  a  non-equilibrium  type.  It  is  interesting  to  point  out 
that  despite  the  fact  that  no  viscous-dissipation  term  was 
added  to  the  Newton’s  equations  of  motion,  the  results 
displayed  in  Fig.  4a,  b  show  some  of  the  defining  features  of 
shock- waves  when  they  are  analyzed  in  the  continuum-level 
simulations  (in  the  presence  of  viscous  dissipation).  Spe¬ 
cifically,  a  steady  shock  is  generated  and  the  shock  width 
decreases  with  an  increase  in  the  shock  strength.  This  point 
has  been  demonstrated  in  Fig.  4d  in  which  curves  labeled 
3,4,  and  5  from  Fig.  4a,  b  have  been  re-plotted. 

The  examination  of  the  results  displayed  in  Figs.  4a-d 
further  reveals  that  in  the  weak-to-intermediate  shock 
strength  range  examined,  the  observed  shock- width  is  quite 
small  (ca.  2-3  nm).  This  finding  suggests  that  the  dominant 
energy  dissipation  processes  captured  by  the  molecular- 
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level  simulations  are  quite  weak.  As  will  be  shown  later, 
these  processes  involve  changes  in  the  soda-lime  glass 
atomic  coordination  and  in  the  size  of  the  smallest  Si-0 
rings.  In  other  words,  as  pointed  out  by  one  of  the 
reviewers  of  this  article,  viscous  dissipation  processes  with 
characteristic  times  at  the  microsecond  time  scale  are  not 
active  in  the  present  simulations.  This,  in  general,  only 
affects  the  width  of  the  shock  profile  but  not  the  values  of 
the  shocked-material  state  variables. 

Determination  of  shock-Hugoniot  relations 


the  present  computational  work.  These  disagreements 
reveal  that  shortcomings/limitations  of  the  present  com¬ 
putational  approach. 

The  Us  vs.  uv  relation  mentioned  above  is  a  simple 
projection  of  the  Hugoniot  to  the  Us-up  plane.  In  the  case 
of  planar  shocks,  of  interest  in  this  study,  the  other  com¬ 
monly  used  Hugoniot  relations  include:  axial  stress,  £n,  vs. 
density,  p  (or  specific  volume,  v  =  1/p);  (mass-based) 
internal  energy  density,  e ,  vs.  p  (or  v);  tn  vs.  uv  and 
temperature  T  vs.  p  (or  v).  These  relations  were  determined 
in  this  study  using  two  distinct  methods: 


The  results  presented  in  Fig.  4a,  b  reveal  the  steady- shock 
profile  and  can  be  used  to  find  a  functional  relationship 
between  the  shock  speed,  Us,  and  the  particle  velocity,  uv. 
This  was  done  in  this  study  and  the  results  obtained 
(depicted  using  discrete  symbols)  along  with  their  least- 
square  linear  fit  (depicted  as  a  straight  line)  are  displayed  in 
Fig.  5.  Also  displayed  in  Fig.  5  is  the  corresponding 
experimental  data  from  [22].  Examination  of  Fig.  5  reveals 
that  the  agreement  between  the  present  computational  data 
and  their  experimental  counterparts  is  fair.  In  particular, 
there  are  two  main  points  of  disagreement  between  the  two 
sets  of  results:  (a)  the  experimental  data  in  the  low  particle 
velocity  (i.e.,  at  low  stress)  region  show  an  anomalous 
material  behavior  which  is  characterized  by  negative  par¬ 
ticle-velocity  dependency  of  the  shock  (or  more  precisely  a 
smooth  wave)  speed;  and  (b)  in  the  high  particle- velocity 
range,  the  experimental  results  show  a  normal  material 
behavior  characterized  by  a  positive  particle-rate  depen¬ 
dency  of  the  shock  speed.  However,  the  magnitude  of  this 
dependency  is  somewhat  larger  than  the  one  predicted  by 
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Fig.  5  The  Hugoniot  relation  pertaining  to  the  particle  velocity 
dependence  of  the  shock  speed  (Color  figure  online) 
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(a)  The  first  method  is  based  on  the  three  so-called  jump 
equations  which  are  defined  as  [30]: 

p~Us  =  p+(Us-up)  (2) 

tb+p~U*  =  t+  +p+(us-up)2  (3) 

+  —  +  O.5E/3  =  e+  +  +  0.5  (t/s  —  up)2  (4) 

These  equations  relate  the  known  downstream  material 
states  (denoted  by  a  superscript  “  —  ”)  and  the  unknown 
upstream  material  states  (denoted  by  a  superscript  “+”) 
associated  with  the  shock  of  a  given  strength  (as  quantified 
by  the  shock  speed  or  the  downstream-to-upstream  particle 
velocity  jump).  These  equations  are  next  combined  with  the 
previously  determined  Us  vs.  up  relation  and  the  prescribed 
(shock-strength  defining  quantity)  Us  or  uv  to  solve  for  the 
unknown  upstream  material  states.  It  should  be  noted  that 
this  method  enables  determination  of  only  material 
mechanical  state  variables  (tn,  e ,  v(=l/p),  Us  and  uv).  To 
obtain  temperature,  a  separate  set  of  equilibrium  NVE  (E — 
total  energy  of  the  system)  molecular  dynamics  simulations 
is  carried  out.  In  each  case,  a  local  computational  sub-cell  is 
defined  containing  only  the  upstream  (shocked)  material. 
The  number  of  particles,  the  volume  of  the  sub-cell  and  its 
total  internal  energy  are  all  maintained  constant.  The 
associated  “equilibrium”  temperature  is  then  calculated 
using  the  time-averages  of  the  atomic  velocities  (see  Eq.  6 
below);  and 

(b)  Time  averages  of  the  atomic  positions,  rb  velocities,  vh 
and  interaction  forces,  ft  (i  is  the  atomic  label)  are  used 
to  compute  the  unknown,  local  (bin-based)  thermo¬ 
mechanical  quantities  using  the  following  standard 
statistical  thermodynamic  relations: 


1 


avg 


(5) 
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T  = 


y3iVbin^ 


-^bin 

i=  1 


avg 


^  /  -^bin 

fn  =  NbinkbT  +  ^  r;  0/;- 
VWn  V  ,=i 


avg 


(6) 

(7) 

(8) 


where  subscript  avg  denotes  time  averaging,  Afbin  and  Vbin 
are  the  number  of  atoms  within  and  the  volume  of  the  bin, 
respectively,  kb  is  the  Boltzmann’s  constant,  ifxotai  is  the 
potential  energy  of  a  system  in  a  given  configuration,  while 
and  “(g)”  indicate  dot  product  and  tensorial  product 
operators,  respectively.  It  should  be  noted  that  the  stress 
relation,  Eq  7  is  a  particular  representation  of  the  standard 
virial-based  stress  equation  [31]. 


It  should  be  noted  that,  in  this  case,  the  bins  were 
defined  within  a  reference  frame  which  is  attached  to  (and 
moves  with)  the  steady-shock  front.  Clearly,  in  this  case 
different  atoms  may  reside  within  a  given  bin  at  different 


simulation  times.  On  the  other  hand,  the  bins  correctly 
collect  the  information  about  the  atoms  (temporarily) 
residing  in  a  given  segment  of  the  steady- shock  profile.  In 
other  words,  time  averages  are  calculated  not  for  a  fixed 
assembly  of  atoms,  but  rather  for  a  transient  set  of  atoms 
associated  with  a  given  moving  bin.  It  should  be  also  noted 
that  since  one  of  the  main  objectives  of  this  study  was 
determination  of  the  Hugoniot  relations,  only  the  data 
pertaining  to  the  bins  located  in  the  fully  shocked  upstream 
region  are  collected  and  analyzed  (for  different  shock- 
strength  conditions).  In  other  words,  the  data  collected  by 
the  bins  located  within  the  shock  profile  and  downstream  of 
the  shock  are  ignored. 

The  results  of  these  two  procedures  are  displayed  in 
Fig.  6a-d.  In  each  of  these  figures,  two  cases  are  shown 
and  labeled  as  “Method  (a)”  or  “Method  (b)”  to  indicate 
the  method  used  for  generation  of  the  corresponding 
results.  The  results  displayed  in  these  figures  show  that,  in 
the  weak-shock  regime,  method  (a)  consistently  over-pre- 
dicts  (by  as  much  as  100%)  the  values  of  the  material  state 
variables  relative  to  method  (b).  This  relative  difference 
continuously  diminishes  with  shock  strength  so  that  in  the 


Fig.  6  a  Axial  stress  versus  specific  volume;  b  energy  versus  specific  volume;  c  axial  stress  versus  particle  velocity  Hugoniot  relations.  Please 
see  text  for  explanation  of  the  “Method  (a)”  and  “Method  (b)”  (Color  figure  online) 
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strong- shock  regime  the  predictions  of  the  two  methods  are 
different  by  less  than  10%. 

The  Hugoniot  relations  displayed  in  Figs.  5  and  6a-d  are 
typically  used  within  a  continuum-level  computational 
analysis  of  shock-wave  generation/propagation  in  two 
ways: 

(a)  They  are  directly  used  in  the  analysis  of  shock-wave 
propagation  under  uniaxial  strain  conditions  [32,  33]; 
and 

(b)  Alternatively,  they  can  be  used  to  derive  a  continuum- 
level  material  model  which  is  consistent  with  the 
material  mechanical  response  under  high-rate,  large- 
strain,  high-pressure  conditions.  Such  a  model  is 
subsequently  used  in  general  three-dimensional,  non¬ 
linear  dynamics  computational  analyses  [22]. 

In  our  recent  work  [6],  it  was  shown  that  the  Hugoniot 
relations  can  be  generated  by  converting  the  corresponding 
isotherms  (obtained  via  quasi- static,  molecular-level 
mechanical  tests).  This  procedure  was  found  to  be  associ¬ 
ated  with  a  number  of  challenges  (e.g.,  a  particular  form  of 
the  equation  of  state  had  to  be  assumed,  several  material 
properties/relations  had  to  be  assessed  independently,  etc.). 
Most  of  these  challenges  were  not  encountered  in  this  study 
since  the  Hugoniot  relations  are  derived  more  directly  from 
the  molecular-level  computational  results. 

Shock-induced  material- state  changes 

The  results  presented  and  discussed  in  the  previous  sections 
clearly  revealed  the  formation  and  propagation  of  planar 
shocks  in  soda-lime  glass  and  enabled  formulation  of  the 
appropriate  shock-Hugoniot  relations.  In  the  present  sec¬ 
tion,  a  more  detailed  investigation  is  carried  out  of  the 
molecular-level  material  microstructure  in  the  wake  of  a 
propagating  planar  shock. 

It  should  be  first  realized  that  the  analysis  of  material 
microstructure  and  its  evolution  due  to  shock  loading  in 
soda-lime  glass  is  quite  challenging  due  to  the  absence  of  a 
crystal  structure  in  this  material.  Namely,  when  molecular- 
level  simulations  of  shock  generation/propagation  are  car¬ 
ried  out  in  (nearly  perfect)  single  crystal  solids  [17,  20],  the 
presence  of  a  crystal  lattice  greatly  facilitates  the  investi¬ 
gation  of  deviations  from  the  long  range  order  (i.e.,  for¬ 
mation  of  various  point,  line  and  planar  defects)  and  the 
nature  of  the  associated  inelastic  deformation  processes. 
Soda-lime  glass  is,  on  the  other  hand,  an  amorphous 
material  in  its  initial  condition  and  remains  so  after  being 
subjected  to  shock  loading.  To  address  the  challenge  of 
material  microstructure  characterization  and  its  changes 
resulting  from  shock  loading,  the  following  microstructural 
parameters  were  monitored:  (a)  the  random  network 
X  parameter.  The  R  parameter  was  not  monitored  since  it  is 


microstructure-insensitive  while  Y  parameter  was  not 
monitored  since  X  +  Y  =  4.0  in  pure  silicate  glasses  (like 
soda  lime  glass);  (b)  the  Si-atom  average  coordination 
number(i.e.,  bonding  structure);  (c)  the  size  of  the  smallest 
Si-0  ring;  and  (d)  the  material’s  average  density. 

In  the  case  of  single-crystalline  solids,  previous  molec¬ 
ular-level  shock- wave  formation/calculation  work  (e.g., 
[17,  20])  established  that  the  steady- wave  condition  is 
attained  not  as  a  result  of  (velocity-dependent)  viscous 
dissipation  (as  is  the  case  for  shocks  in  fluids),  but  rather  as 
a  result  of  inelastic  deformation  (permanent  slippage  of 
crystal  planes  and  the  formation  of  crystal  defects)  and 
phase  transformation  processes.  The  microstructural 
parameters  identified  above  are  used  to  characterize  the 
type  and  the  extent  of  these  processes. 

X  parameter  evaluation 

In  general,  changes  in  the  X  parameter  are  the  result  of 
competition  between  Si  and  O  bond  breaking  processes 
(increases  the  average  number  of  non-bridging  oxygens 
per  un-shocked  material  polyhedron)  and  the  processes 
which  lead  to  an  increase  in  Si-atom  coordination  number 
(decrease  the  X  parameter  value).  Since  the  instances  of 
Si-atom  coordination-number  increase  are  found  to  greatly 
out-number  the  Si-0  bond-breaking  instances,  the  value  of 
the  X  parameter  has  been  found  to  decrease  by  5-10%  from 
its  equilibrium  value  (=0.81). 

Si-atom  average  coordination  number 

Examination  of  the  material  micro  structure  in  the  shock 
wake  revealed  that  the  Si-atom  average  coordination 
number  increases  from  its  initial  value  of  4.0  by  2-3%.  An 
example  of  the  shock-loading  induced  increase  in  Si-atom 
coordination  number  is  depicted  in  Fig.  7a-d.  To  aid  in 
visualization/interpretation  of  the  microstructural/topolo- 
gical  changes  experienced  by  soda-lime  glass  during  shock 
loading,  only  a  20-30  atom  exemplary  region  of  the 
upstream  material  is  analyzed  in  these  figures.  In  addition, 
sodium  cations  are  not  displayed. 

In  Fig.  7a,  c,  x-y  and  y-z  projections  are  given  of  the 
material  region  under  consideration  in  its  initial  state.  The 
corresponding  projections  for  the  same  material  region  after 
the  passing  of  a  shock  are  displayed  in  Fig.  7b,  d.  For 
improved  clarity,  three  silicon  atoms  are  highlighted  in 
green.  To  make  apparent  the  Si-atom  coordination-number 
increase,  one  of  the  Si  atoms  is  labeled  as  “57”  along  with  the 
oxygen  atoms  bonded  to  it.  It  is  clear  that  the  silicon  atom  in 
question  was  initially  four-fold  coordinated  and  it  became 
five-fold  coordinated  upon  the  passage  of  the  shock. 

It  should  be  noted  that  the  molecular-level  microstruc¬ 
tural  changes  described  above  are  a  manifestation  of  the 
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Fig.  7  Development  of  five¬ 
fold  coordinated  silicon  atoms 
in  soda-lime  glass  under  shock 
loading:  a  material  initial  state, 
x-y  projection;  b  post- shock 
state  of  the  same  material  region 
as  in  (a);  c  and  d  y-z  projection 
of  (a)  and  (b),  respectively 
(Color  figure  online) 


transition  of  soda-lime  glass  to  a  state  that  is  energetically 
favored  at  high  shock-induced  stresses  (large  densities). 
This  finding  is  consistent  with  the  experimental  observa¬ 
tion  of  Alexander  et  al.  [22]  who  reported  the  formation  of 
stishovite-like  domains  containing  six-fold  coordinated 
silicon  atoms  in  soda-lime  glass  at  ca.  30  GPa  shock- 
induced  axial  stress  levels. 

Si-0  ring  minimum  size 

It  should  be  first  noted  that  a  comprehensive  statistics- 
based  Si-0  ring  minimum-size  analysis  is  quite  compli¬ 
cated  and  is  beyond  the  scope  of  this  study.  Such  an 
analysis  will  be  carried  out  in  our  future  work.  Instead, 
only  a  qualitative  analysis  of  the  molecular-level  shocked 
material  microstructure  is  carried  out  here  to  obtain  evi¬ 
dence  for  the  potential  reduction  in  the  minimum  size  of 
the  Si-0  rings.  An  example  of  two  Si-0  rings  each  con¬ 
taining  five  silicon  atoms  found  in  the  as-shocked  material 
is  shown  in  Fig.  8.  For  improved  clarity,  these  rings  are 
highlighted  in  purple  and  yellow.  Similar  rings  were  not 


observed  in  the  initial  un-shocked  material.  This  finding 
was  reconfirmed  using  different  (yet,  all  equilibrium)  initial 
molecular  structures.  Thus,  this  study  suggests  that  shock 
loading  causes  a  reduction  in  the  minimum  Si-0  ring  size. 

Material’s  average  density 

It  should  be  also  recognized  that  shock  loading  leads  to  a 
permanent  2-4%  (shock  strength-dependent)  increase  in 
the  material  density.  This  material  densification  process 
along  with  the  aforementioned  processes  leading  to  an 
increase  in  the  silicon-atom  coordination  and  smaller  Si-0 
ring  formation  are  all  associated  with  energy  absorption/ 
dissipation  and,  hence,  are  expected  to  play  an  important 
role  in  the  blast/ballistic  impact  mitigation  potential  of 
soda-lime  glass. 

Final  remarks 

Within  this  study,  molecular-level  computational  methods 
are  employed  to  study  various  phenomena  accompanying 
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Fig.  8  Observations  of  smaller  Si-O  rings  (containing  five  Si  atoms) 
in  the  as-shocked  soda-lime  glass.  Similar  small-size  rings  Si-O  rings 
were  not  observed  in  the  initial/un- shocked  material  state.  The  rings 
in  question  have  been  highlighted  for  clarity  (Color  figure  online) 

shock- wave  generation  and  propagation  in  soda-lime  glass. 
It  should  be  noted  that  this  study  does  not  suggest  that 
molecular-level  analyses  of  shock  generation  and  propa¬ 
gation  should  replace  the  corresponding  continuum-level 
(hydrodynamic)  analyses.  The  latter  are  far  better  suited 
(and  feasible)  for  studying  the  behavior  of  real-life  engi¬ 
neering  systems  (e.g.,  vehicle  windshield  subjected  to  blast 
impact).  Rather,  the  present  approach  is  highly  beneficial 
relative  to  the  identification  and  characterization  of  the 
nanoscale  phenomena/processes  taking  place  at  the  shock 
front.  Once  these  phenomena/processes  are  well  under¬ 
stood  and  characterized  (a  formidable  task)  the  knowledge 
gained  can  be  used  to  formulate  (and  parameterize)  more 
physically  based  material  models  suitable  for  use  in  con¬ 
tinuum-level  analyses. 

Summary  and  conclusions 

Based  on  the  results  obtained  in  this  study,  the  following 
summary  remarks  and  main  conclusions  can  be  drawn: 

1.  Various  phenomena  accompanying  the  formation  and 
propagation  of  a  planar  shock-wave  within  soda-lime 
glass,  a  material  commonly  used  in  transparent  armor 
applications,  are  investigated  using  molecular-level 
computational  methods. 

2.  The  results  obtained  show  that  even  without  the  use  of 
a  viscous-dissipation-based  thermostat,  a  steady-wave 
planar  shock  profile  can  readily  be  established  in  this 
material. 

3.  The  time-averaged  results  pertaining  to  the  atomic 
positions,  velocities  and  interaction  forces  are  used  to 


construct  the  appropriate  shock-Hugoniot  relations, 
the  relations  which  define  the  locus  of  stress,  energy 
density,  mass  density,  temperature,  and  particle  veloc¬ 
ity  in  soda-lime  glass  swept  by  a  shock  propagating  at 
a  given  speed. 

4.  Detailed  examination  of  the  molecular-level  micro¬ 
structure  evolution  in  the  shock- wave  wake  is  carried 
out  to  identify  the  nature  of  energy-absorbing  and  shock- 
wave  spreading  mechanisms.  The  results  revealed  that 
shock  loading  causes  extensive  changes  in  atomic 
coordination  and  the  bond  structure  as  well  as  a  2-4% 
(shock  strength-dependent)  density  increase.  Specifi¬ 
cally,  the  atomic  coordination  of  many  silicon  atoms  has 
been  found  to  increase  from  four  to  five  and  numerous 
smaller  Si-0  rings  are  observed.  These  processes 
are  associated  with  substantial  energy  absorption  and 
dissipation  and  are  believed  to  greatly  influence 
the  blast/ballistic  impact  mitigation  potential  of  soda- 
lime  glass. 
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